home *** CD-ROM | disk | FTP | other *** search
/ QRZ! Ham Radio 8 / QRZ Ham Radio Callsign Database - Volume 8.iso / pc / files / ant_nec / nec81tar.z / nec81tar / factr.f < prev    next >
Text File  |  1991-05-13  |  8KB  |  354 lines
  1. C $TITLE: 'FACTR'
  2. C $NOFLOATCALLS
  3. C
  4. C
  5.       SUBROUTINE FACTR(A,D,N,NDIM,IP,LD2)
  6. C
  7. C     SUBROUTINE TO FACTOR A MATRIX INTO A UNIT LOWER TRIANGULAR MATRIX
  8. C     AND AN UPPER TRIANGULAR MATRIX USING THE GAUSS-DOOLITTLE ALGORITHM
  9. C     PRESENTED ON PAGES 411-416 OF A. RALSTON--A FIRST COURSE IN
  10. C     NUMERICAL ANALYSIS.  COMMENTS BELOW REFER TO COMMENTS IN RALSTONS
  11. C     TEXT.    (MATRIX TRANSPOSED.
  12. C
  13.       INTEGER*4 N,NDIM
  14.       INTEGER R,RM1,RP1,PJ,PR
  15.       REAL*8 ELMAG,DMAX
  16. CLARGE: A
  17.       COMPLEX A
  18.       COMPLEX*16 D,ARJ,SUM
  19.       DIMENSION D(LD2)
  20. C**
  21. C**      DIMENSION IP(NDIM),A(NDIM,NDIM)
  22.       DIMENSION IP(NDIM),A(NDIM,1)
  23. C**
  24. C     D      WRITE(*,*) '   FACTR: START N=',N,' NDIM=',NDIM
  25. C**
  26.       IFLG=0
  27. C**
  28.       SUM=DCMPLX(0.,0.)
  29.       DO 9 R=1,N
  30. C
  31. C     STEP 1
  32. C
  33.       DO 1 K=1,N
  34.       D(K)=A(R,K)
  35. 1     CONTINUE
  36. C
  37. C     STEPS 2 AND 3
  38. C
  39.       RM1=R-1
  40.       IF (RM1.LT.1) GO TO 4
  41.       DO 3 J=1,RM1
  42.       PJ=IP(J)
  43.       ARJ=D(PJ)
  44.       A(R,J)=ARJ
  45.       D(PJ)=D(J)
  46.       JP1=J+1
  47.       DO 2 I=JP1,N
  48.       D(I)=D(I)-A(J,I)*ARJ
  49. 2     CONTINUE
  50. 3     CONTINUE
  51. 4     CONTINUE
  52. C
  53. C     STEP 4
  54. C
  55.       DMAX=DREAL(D(R)*DCONJG(D(R)))
  56.       IP(R)=R
  57.       RP1=R+1
  58.       IF (RP1.GT.N) GO TO 6
  59.       DO 5 I=RP1,N
  60.       ELMAG=DREAL(D(I)*DCONJG(D(I)))
  61.       IF (ELMAG.LT.DMAX) GO TO 5
  62.       DMAX=ELMAG
  63.       IP(R)=I
  64. 5     CONTINUE
  65. 6     CONTINUE
  66.       IF (DMAX.LT.1.D-10) IFLG=1
  67.       PR=IP(R)
  68.       A(R,R)=D(PR)
  69.       D(PR)=D(R)
  70. C**
  71.       SUM=SUM+A(R,R)
  72. C**
  73. C
  74. C     STEP 5
  75. C
  76.       IF (RP1.GT.N) GO TO 8
  77.       ARJ=1.D0/A(R,R)
  78.       DO 7 I=RP1,N
  79.       A(R,I)=D(I)*ARJ
  80. 7     CONTINUE
  81. 8     CONTINUE
  82.       IF (IFLG.EQ.0) GO TO 9
  83.       WRITE(*,10)  R,DMAX
  84.       IFLG=0
  85. 9     CONTINUE
  86.       RETURN
  87. C
  88. 10    FORMAT (' PIVOT(',I3,')=',1P,D16.8)
  89.       END
  90. C
  91. C
  92. C
  93.       SUBROUTINE SOLVES(A,B,Y,NEQ,NP,N,MP,M,IP,NRH,IFL1,IFL2,LD2,
  94.      1 IRESRV)
  95. C
  96. C     SUBROUTINE SOLVES, FOR SYMMETRIC STRUCTURES, HANDLES THE
  97. C     TRANSFORMATION OF THE RIGHT HAND SIDE VECTOR AND SOLUTION OF THE
  98. C     MATRIX EQ.
  99. C
  100.       COMMON/SMAT/ SSX(16,16)
  101.       INTEGER*4 NPEQ,NROW,NEQ,NP,N,MP,M
  102.       INTEGER*4 IMAT,NPBLK,NLAST,NLSYM,NBBX,NPBX,NLBX,NBBL,NPBL,NLBL
  103.       REAL*8 FNORM
  104. CLARGE: A,B
  105.       COMPLEX A,B
  106.       COMPLEX*16 Y,SUM
  107.       COMPLEX*16 SSX
  108.       DIMENSION A(1),IP(1),B(NEQ,NRH),Y(LD2)
  109.       COMMON/MATPAR/ICASE,NBLOKS,NPBLK,NLAST,NBLSYM,NPSYM,NLSYM,IMAT,
  110.      1 ICASX,NBBX,NPBX,NLBX,NBBL,NPBL,NLBL
  111. C**
  112. C     D      WRITE(*,*) '   SOLVES: START IFL1=',IFL1,' IFL2=',IFL2
  113.       IRESRV=IRESRV
  114. C**
  115.       NPEQ=NP+2*MP
  116.       NOP=NEQ/NPEQ
  117.       FNOP=NOP
  118.       FNORM=1.D0/FNOP
  119.       NROW=NEQ
  120.       IF (ICASE.GT.3) NROW=NPEQ
  121.       IF (NOP.EQ.1) GO TO 11
  122.       DO 10 IC=1,NRH
  123.       IF (N.EQ.0.OR.M.EQ.0) GO TO 6
  124.       DO 1 I=1,NEQ
  125. 1     Y(I)=B(I,IC)
  126.       KK=2*MP
  127.       IA=NP
  128.       IB=N
  129.       J=NP
  130.       DO 5 K=1,NOP
  131.       IF (K.EQ.1) GO TO 3
  132.       DO 2 I=1,NP
  133.       IA=IA+1
  134.       J=J+1
  135. 2     B(J,IC)=Y(IA)
  136.       IF (K.EQ.NOP) GO TO 5
  137. 3     DO 4 I=1,KK
  138.       IB=IB+1
  139.       J=J+1
  140. 4     B(J,IC)=Y(IB)
  141. 5     CONTINUE
  142. C
  143. C     TRANSFORM MATRIX EQ. RHS VECTOR ACCORDING TO SYMMETRY MODES
  144. C
  145. 6     DO 10 I=1,NPEQ
  146.       DO 7 K=1,NOP
  147.       IA=I+(K-1)*NPEQ
  148. 7     Y(K)=B(IA,IC)
  149.       SUM=Y(1)
  150.       DO 8 K=2,NOP
  151. 8     SUM=SUM+Y(K)
  152.       B(I,IC)=SUM*FNORM
  153.       DO 10 K=2,NOP
  154.       IA=I+(K-1)*NPEQ
  155.       SUM=Y(1)
  156.       DO 9 J=2,NOP
  157. 9     SUM=SUM+Y(J)*DCONJG(SSX(K,J))
  158. 10    B(IA,IC)=SUM*FNORM
  159. 11    IF (ICASE.LT.3) GO TO 12
  160.       REWIND IFL1
  161.       REWIND IFL2
  162. C
  163. C     SOLVE EACH MODE EQUATION
  164. C
  165. 12    DO 16 KK=1,NOP
  166.       IA=(KK-1)*NPEQ+1
  167.       IB=IA
  168.       IF (ICASE.NE.4) GO TO 13
  169.       I=NPEQ*NPEQ
  170.       READ (IFL1) (A(J),J=1,I)
  171.       IB=1
  172. 13    IF (ICASE.EQ.3.OR.ICASE.EQ.5) GO TO 15
  173.       DO 14 IC=1,NRH
  174. 14    CALL SOLVE(A(IB),B(IA,IC),Y,NPEQ,NROW,IP(IA),LD2)
  175.       GO TO 16
  176. 15    CALL LTSOLV(A,B(IA,1),Y,IP(IA),NPEQ,NEQ,NRH,IFL1,IFL2,LD2)
  177. 16    CONTINUE
  178. C**
  179. C     D     IF (NOP.EQ.1) WRITE(*,*) '   SOLVES: EARLY RETURN'
  180. C**
  181.       IF (NOP.EQ.1) RETURN
  182. C
  183. C     INVERSE TRANSFORM THE MODE SOLUTIONS
  184. C
  185.       DO 26 IC=1,NRH
  186.       DO 20 I=1,NPEQ
  187.       DO 17 K=1,NOP
  188.       IA=I+(K-1)*NPEQ
  189. 17    Y(K)=B(IA,IC)
  190.       SUM=Y(1)
  191.       DO 18 K=2,NOP
  192. 18    SUM=SUM+Y(K)
  193.       B(I,IC)=SUM
  194.       DO 20 K=2,NOP
  195.       IA=I+(K-1)*NPEQ
  196.       SUM=Y(1)
  197.       DO 19 J=2,NOP
  198. 19    SUM=SUM+Y(J)*SSX(K,J)
  199. 20    B(IA,IC)=SUM
  200.       IF (N.EQ.0.OR.M.EQ.0) GO TO 26
  201.       DO 21 I=1,NEQ
  202. 21    Y(I)=B(I,IC)
  203.       KK=2*MP
  204.       IA=NP
  205.       IB=N
  206.       J=NP
  207.       DO 25 K=1,NOP
  208.       IF (K.EQ.1) GO TO 23
  209.       DO 22 I=1,NP
  210.       IA=IA+1
  211.       J=J+1
  212. 22    B(IA,IC)=Y(J)
  213.       IF (K.EQ.NOP) GO TO 25
  214. 23    DO 24 I=1,KK
  215.       IB=IB+1
  216.       J=J+1
  217. 24    B(IB,IC)=Y(J)
  218. 25    CONTINUE
  219. 26    CONTINUE
  220. C**
  221. C     D      WRITE(*,*) '   SOLVES: RETURN AT END'
  222. C**
  223.       RETURN
  224.       END
  225. C
  226. C
  227. C
  228.       SUBROUTINE SOLVE(A,B,Y,N,NDIM,IP,LD2)
  229. C
  230. C     SUBROUTINE TO SOLVE THE MATRIX EQUATION LU*X=B WHERE L IS A UNIT
  231. C     LOWER TRIANGULAR MATRIX AND U IS AN UPPER TRIANGULAR MATRIX BOTH
  232. C     OF WHICH ARE STORED IN A.  THE RHS VECTOR B IS INPUT AND THE
  233. C     SOLUTION IS RETURNED THROUGH VECTOR B.    (MATRIX TRANSPOSED.
  234. C
  235.       INTEGER*4 N,NDIM
  236.       INTEGER PI
  237. CLARGE: A,B
  238.       COMPLEX A,B
  239.       COMPLEX*16 Y,SUM
  240.       DIMENSION A(NDIM,1),B(NDIM),IP(NDIM),Y(LD2)
  241. C**
  242. C     D      WRITE(*,*) '   SOLVE: START'
  243. C**
  244. C
  245. C     FORWARD SUBSTITUTION
  246. C
  247.       DO 3 I=1,N
  248.       PI=IP(I)
  249.       Y(I)=B(PI)
  250.       B(PI)=B(I)
  251.       IP1=I+1
  252.       IF (IP1.GT.N) GO TO 2
  253.       DO 1 J=IP1,N
  254.       B(J)=B(J)-A(I,J)*Y(I)
  255. 1     CONTINUE
  256. 2     CONTINUE
  257. 3     CONTINUE
  258. C
  259. C     BACKWARD SUBSTITUTION
  260. C
  261.       DO 6 K=1,N
  262.       I=N-K+1
  263.       SUM=DCMPLX(0.,0.)
  264.       IP1=I+1
  265.       IF (IP1.GT.N) GO TO 5
  266.       DO 4 J=IP1,N
  267.       SUM=SUM+A(J,I)*B(J)
  268. 4     CONTINUE
  269. 5     CONTINUE
  270.       B(I)=(Y(I)-SUM)/A(I,I)
  271. 6     CONTINUE
  272. C**
  273. C     D      WRITE(*,*) '   SOLVE: RETURN'
  274. C**
  275.       RETURN
  276.       END
  277. C
  278.       SUBROUTINE LTSOLV (A,B,Y,IX,NROW,NEQ,NRH,IFL1,IFL2,LD2)
  279. C
  280. C     LTSOLV SOLVES THE MATRIX EQ. Y(R)*LU(T)=B(R) WHERE (R) DENOTES ROW
  281. C     VECTOR AND LU(T) DENOTES THE LU DECOMPOSITION OF THE TRANSPOSE OF
  282. C     THE ORIGINAL COEFFICIENT MATRIX.  THE LU(T) DECOMPOSITION IS
  283. C     STORED ON TAPE 5 IN BLOCKS IN ASCENDING ORDER AND ON FILE 3 IN
  284. C     BLOCKS OF DESCENDING ORDER.
  285. C
  286.       INTEGER*4 I2,I,J,K,NROW,NEQ,IDM1
  287.       INTEGER*4 IMAT,NPBLK,NLAST,NLSYM,NBBX,NPBX,NLBX,NBBL,NPBL,NLBL
  288. CLARGE A,B
  289.       COMPLEX A,B
  290.       COMPLEX*16 Y,SUM
  291.       COMMON/MATPAR/ICASE,NBLOKS,NPBLK,NLAST,NBLSYM,NPSYM,NLSYM,IMAT,
  292.      1 ICASX,NBBX,NPBX,NLBX,NBBL,NPBL,NLBL
  293.       DIMENSION A(NROW,1),B(NEQ,NRH),Y(LD2),IX(NEQ)
  294. C
  295. C     FORWARD SUBSTITUTION
  296. C
  297. C     D      WRITE(*,*) '    LTSOLV: START'
  298.       IDM1=1
  299.       I2=2*NPSYM*NROW
  300.       DO 4 IXBLK1=1,NBLSYM
  301.       CALL BLCKIN (A,IDM1,I2,1,121,IFL1)
  302.       K2=NPSYM
  303.       IF (IXBLK1.EQ.NBLSYM) K2=NLSYM
  304.       JST=(IXBLK1-1)*NPSYM
  305.       DO 4 IC=1,NRH
  306.       J=JST
  307.       DO 3 K=1,K2
  308.       JM1=J
  309.       J=J+1
  310.       SUM=(0.,0.)
  311.       IF (JM1.LT.1) GO TO 2
  312.       DO 1 I=1,JM1
  313. 1     SUM=SUM+A(I,K)*B(I,IC)
  314. 2     B(J,IC)=(B(J,IC)-SUM)/A(J,K)
  315. 3     CONTINUE
  316. 4     CONTINUE
  317. C
  318. C     BACKWARD SUBSTITUTION
  319. C
  320.       JST=NROW+1
  321.       DO 8 IXBLK1=1,NBLSYM
  322.       CALL BLCKIN (A,IDM1,I2,1,122,IFL2)
  323.       K2=NPSYM
  324.       IF (IXBLK1.EQ.1) K2=NLSYM
  325.       DO 7 IC=1,NRH
  326.       KP=K2+1
  327.       J=JST
  328.       DO 6 K=1,K2
  329.       KP=KP-1
  330.       JP1=J
  331.       J=J-1
  332.       SUM=DCMPLX(0.,0.)
  333.       IF (NROW.LT.JP1) GO TO 6
  334.       DO 5 I=JP1,NROW
  335. 5     SUM=SUM+A(I,KP)*B(I,IC)
  336.       B(J,IC)=B(J,IC)-SUM
  337. 6     CONTINUE
  338. 7     CONTINUE
  339. 8     JST=JST-K2
  340. C
  341. C     UNSCRAMBLE SOLUTION
  342. C
  343.       DO 10 IC=1,NRH
  344.       DO 9 I=1,NROW
  345.       IXI=IX(I)
  346. 9     Y(IXI)=B(I,IC)
  347.       DO 10 I=1,NROW
  348. 10    B(I,IC)=Y(I)
  349. C**
  350. C     D      WRITE(*,*) '    LTSOLV: RETURN'
  351. C**
  352.       RETURN
  353.       END
  354.